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Abstract 

We deal with the Cauchy problem for the space-time fractional diffusion 
equation, which is obtained from the standard diffusion equation by replac- 
ing the second-order space derivative with a Riesz-Feller derivative of order 
a G (0,2] and skewness 9 {\9\ < min{a,2 — a}), and the first-order time 
derivative with a Caputo derivative of order (3 G (0, 2] . The fundamental 
solution (Green function) for the Cauchy problem is investigated with re- 
spect to its scaling and similarity properties, starting from its Fourier-Laplace 
representation. We review the particular cases of space-fractional diffusion 
{0 < a < 2 , /? = 1} , time-fractional diffusion {a = 2 , < /? < 2} , and 
neutral-fractional diffusion {0 < a = /5 < 2} , for which the fundamental 
solution can be interpreted as a spatial probability density function evolving 



in time. Then, by using the Melhn transform, we provide a general repre- 
sentation of the Green functions in terms of Melhn-Barnes integrals in the 
complex plane, which allows us to extend the probability interpretation to the 
ranges {0<a<2}n{0</3<l} and {I < P < a <2}. Furthermore, from 
this representation we derive explicit formulae (convergent series and asymp- 
totic expansions), which enable us to plot the spatial probability densities 
for different values of the relevant parameters a,6, j3 . 

MSG. 26A33, 33E12, 33C40, 44A10, 45K05, 60G52. 

Key Words: diffusion, fractional derivative, Fourier transform, Laplace trans- 
form, Mellin transform, Mittag-Leffler function, Wright function, Mellin- 
Barnes integrals. Green function, stable probability distributions. 

1. Introduction 

A space-time fractional diffusion equation, obtained from the standard diffu- 
sion equation by replacing the second order space-derivative by a fractional 
Riesz derivative order a > and the first order time- derivative by a frac- 
tional derivative of order /5 > (in Caputo or Riemann-Liouville sense), has 
been recently treated by a number of authors, see for example Saichev and 
Zaslavsky [38] , Uchaikin and Zolotarev |18] , Gorenfio, Iskenderov and Luchko 
[T7] . Scalas, Gorenfio and Mainardi [12], Metzler and Klafter [3¥J. For other 
treatments of the space fractional and/or time fractional diffusion equations 
we refer the reader to the references cited therein. See below for the restric- 
tions on a and j3 . In this paper we intend to complement the results obtained 
in [17] by allowing asymmetry in the space fractional derivative. We thus 
consider the space-time fractional diffusion equation 

^D^u{x,t) = tD^u{x,t) , xeR, teM+, (1.1) 

where the a , 6 , (3 are real parameters always restricted as follows 

0<a<2, 1^1 < min{a,2 - a} , < (3 < 2 . (1.2) 

In (1.1) u = u{x,t) is the (real) field variable, ^Dg is the Riesz- Feller space- 
fractional derivative of order a and skewness 6 , and is the Caputo 

time-fractional derivative of order f3 . These fractional derivatives are integro- 
differential operators to be defined later. 
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The paper is divided as follows. In Section 2 we provide the reader 
with the essential notions and notations concerning the Fourier, Laplace and 
Mellin transforms, which are necessary in the following. 

In Section 3 we introduce the Cauchy problem for the equation (1.1) and 
find the corresponding fundamental solution G^^(a;, t) (the Green function) 
in terms of its Fourier- Laplace transform from which we derive its general 
scaling properties and the similarity variable xji^l"- . We shall get the funda- 
mental formula 

Gi^^{x, t) = K^ix/t-^) , 7 = , (1.3) 

where K^^ p{x) is referred to as the reduced Green function. 

In Section 4, we consider the particular cases {0 < a < 2 , /3 = 1} {space 
fractional diffusion), {a = 2 , < /3 < 2} {time fractional diffusion), and 
{0 < q; = /9 < 2} {neutral fractional diffusion), for which the fundamental 
solution can be interpreted as a spatial probability density function {pdf), 
evolving in time. 

In Section 5 we show a composition rule for the Green function, valid for 

{0 < a < 2} n {0 < /3 < 1} 

which ensures its probability interpretation in this range. 

In Section 6, we provide a general representation for the (reduced) Green 
function in terms of a Mellin-Barnes integral in the complex plane, which 
enables us to extend the probability interpretation to the range 

{!< /3 < a < 2}. 

In Section 7, we derive for the Green function explicit formulae (conver- 
gent series and asymptotic expansions), whose form depends on the relation 
between the parameters a and (3 . By means of a suitable matching between 
the convergent and the asymptotic expansion we shall be able to compute 
the Green function in all the cases in which it is interpretable as a probability 
density. 

Finally, Section 8 is devoted to concluding discussions, and a summary 
of the results in which we present plots of the Green function for a number 
of cases. 
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2. Notions and Notations 



For the sake of the reader's convenience here we present an introduction to 
the Riesz-Feller and Caputo fractional derivatives starting from their rep- 
resentation in the Fourier and Laplace transform domain, respectively. So 
doing we avoid the subtleties lying in the inversion of fractional integrals, 
see e.g. |39], [21]. We also recall the main properties of the Mellin transform 
that will be used later. 

Since in what follows we shall meet only real or complex-valued functions 
of a real variable that are defined and continuous in a given open interval 
X = (a, h) , — oo < a < h < +00 , except, possibly, at isolated points where 
these functions can be infinite, we restrict our presentation of the integral 
transforms to the class of functions for which the Riemann improper integral 
on X absolutely converges. In so doing we follow Marichev ^32j and we denote 
this class by L^{I) or L'^(a, h) . 

The Fourier transform and the Riesz-Feller space-fractional derivative 
Let 

^ r+00 . ■ 

f{K)=T{f{x)]K}= e+^'^^/(x)dx, KGiR, (2.1a) 

J —00 

be the Fourier transform of a function /(x) G L'^{1R), and let 

/(x)= J-U/(«:);x} = — / e-'^^f{K)dK, xeM, (2.16) 

Zn J-00 

be the inverse Fourier transfornll]. For a sufficiently well-behaved function 
f{x) we define the Riesz-Feller space-fractional derivative of order a and 
skewness 9 as 

^{.D,"/(x);/t} = -<(«:) /(«:), (2.2) 
^/;^(ft;) = |ft;|°e^(signft:)^vr/2^ 0<a<2, |e| < min {a, 2 - a} . (2.3) 

We note that the allowed region for the parameters a and 6 turns 
out to be a diamond in the plane {a , 6} with vertices in the points 
(0, 0) , (1, 1) , (1, —1) , (2, 0) , that we call the Feller- Takayasu diamond, see 
Fig. 1. Thus, we recognize that the Riesz-Feller derivative is required to be 



^If f(x) is piecewise differentiable, then the formula (2.1b) holds true at all points 
where f{x) is continuous and the integral in it must be understood in the sense of the 
Cauchy principal value. 
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Fig. IThe Feller- Takayasu diamond 



the pseudo-differential operatoiEl whose symbol — ^'^(k) is the logarithm of 
the characteristic function of a general Levy strictly stable probability density 
with index of stability a and asymmetry parameter 6 (improperly called skew- 
ness) according to Feller's parameterization [131 [H] ? as revisited by Gorenfio 
and Mainardi, see [22| [231 [2i] . 

The operator defined by (2.2)-(2.3) has been referred to as the Riesz-Feller 
fractional derivative since it is obtained as the left inverse of a fractional 
integral originally introduced (for 6 = and a; 7^ 1) by Marcel Riesz in the 
late 1940's, known as the Riesz potential, and then generalized (for 6^0) 
by William Feller in 1952, see [ISj, t39j. 

For more details on Levy stable densities we refer the reader to special- 
istic treatises, as Feller [2], Zolotarev [IH], Samorodnitsky and Taqqu [10], 
Janicki and Weron [25], Sato[?l], Uchaikin and Zolotarev [48], where different 
notations are adopted. We like to refer also to the 1986 paper by Schneider 
[33j . where he first provided the Fox if-function representation of the stable 
distributions (with a 7^ 1) and to the 1990 book by Takayasu [H] where he 
first gave the diamond representation in the plane {a, 6}. 

For ^ = we have a symmetric operator with respect to x , which can be 

^Let us recall that a generic pseudo-differential operator A, acting with respect 
to the variable x € JR , is defined through its Fourier representation, namely 
e*"^ A [f{x)] dx = A{k) /(k) , where A{k) is referred to as symbol of A , given as 

A{k) = (ylc-*«^) e+*''"^ . 



interpreted as 



X 



as can be formally deduced by writing — = — (k^)"/^. We thus recog- 
nize that the operator Dq is related to a power of the positive definitive 
operator —^D^ = and must not be confused with a power of the first 
order differential operator r^D = for which the symbol is —in . An alter- 
native illuminating notation for the symmetric fractional derivative is due to 
Zaslavsky, see e.g. [3B], and reads 



In its regularized form valid for < a < 2 the Riesz space-fractional 
derivative admits the explicit represent at iorlf] 

,OSm ^ sin (f ) ^»/(^±i)^|M±A^,,. 

For a = 1 the Riesz derivative is related to the Hilbert transform as first 
noted by Feller in 1952 in his pioneering paper [13], resulting in the formula 



For < a < 2 and |^| < min {a^ 2 — a} the Riesz-Feller derivative reads 

For a = 1 we obtain the composite formula 

.,Dl fix) = [cos(^7r/2) + sin(e7r/2) ^d] f{x) , (2.9) 



■^The representation (2.6), based on a suitable regularization of a hyper-singular inte- 
gral, can be found in Samko, Kilbas & Marichev [39] as formula (12.1') and is more explicit 
and convenient than other ones available in the literature, see e.g. Saichev & Zaslavsky 
[38] , Uchaikin & Zolotarev [48] , in that it is valid in the whole range < a < 2 . Gorenflo 
and Mainardi have used it in [53] , where they have shown that it holds also in the singular 
case a = 1 . 



6 



that in the extremal cases 6 = ±1 reduces to 

^Dl, = ±^D= ±-^. (2.10) 
ax 

The Laplace transform and the Caputo fractional derivative 

Now we present an introduction to the Caputo fractional derivative starting 
from its representation in the Laplace transform domain and contrasting it 
to the standard Riemann-Liouville fractional derivative. 
Let 

fis) = C {fit); s} = e fit) dt , 3fJ (s) > a/ , (2.11a) 

be the Laplace transform of a function fit) G C'^iO, T) , VT > and let 

fit)=C-'{fis);t} = — e'^fis)ds, 3fJ(s)=7>a/, (2.116) 

with t > , be the inverse Laplace transfornfl For a sufficiently well-behaved 
function fit) we define the Caputo time-fractional derivative of order /3 (m — 
l < /5 < m , m G IN) through 

m—l 

C{,DP,fit);s} = s^fis)-J2s^~'-'f^'\0^), m-l<f3<m. (2.12) 

k=0 

This leads to define, see e.g. [S], [2T], 

1 /■* f^'^\T)dT 



Vim - p) Jo it - t)P+^-^ ' 



m — 1 < (3 < m . 



(2.13) 



The operator defined by (2.12)-(2.13) has been referred to as the Caputo 
fractional derivative since it was introduced by Caputo in the late 1960 's for 



sufficient condition of the existence of the Laplace transform is that the original 
function is of exponential order as i oo . This means that some constant a f exists such 
that the product e~°^* 1/(^)1 is bounded for all t greater than some T . Then f{s) exists 
and is analytic in the half plane 5R(s) > a/ . If f{t) is piecewise differentiable, then the 
formula (2.11b) holds true at all points where /(i) is continuous and the (complex) integral 
in it must be understood in the sense of the Cauchy principal value. 
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modelling the energy dissipation in some anelastic materials with memory, 
see in E]- A former review of the theoretical aspects of this derivative with 
applications in viscoelasticity was given in 1971 by Caputo and Mainardi [8], 
with special emphasis to the long-memory effects. 

The reader should observe that the Caputo fractional derivative differs 
from the usual Riemann-Liouville fractional derivative which, defined as the 
left inverse of the Riemann-Liouville fractional integral, is here denoted as 
tD/^ fit) . We have, see e.g. [39], 



dt"" 
df^ 



fjr) dr 
T{m-p) Jo (t-r)^+i-'^ 



m — 1 < j3 < m, 



(3 = m . 



(2.14) 



When the order is not integer, Gorenfio and Mainardi have shown the follow- 
ing relationships between the two fractional derivatives (when both of them 
exist), see e.g. [21] . 



tD^. fit) = tD^ 



m—l 



/W-E/^'Ho+)^ 

fc=0 



m — 1 < [3 < m . 



(2.15) 



or 



m— 1 



m-1 < P <m. (2.16) 



The Caputo fractional derivative is of course more restrictive than the 
Riemann-Liouville fractional derivative in that the derivative of order m is 
required to exist and be absolutely Laplace transformable. Whenever we use 
the operator , we (tacitly) assume that this condition is met. 

The Caputo fractional derivativ^, practically ignored in the mathemat- 
ical treatises, represents a sort of regularization in the time origin for the 
Riemann-Liouville fractional derivative. Recently, it has been extensively 
investigated by Gorenfio & Mainardi [25 and by Podlubny [SB] in view of 
its major utility in treating physical and engineering problems which require 
standard initial conditions. Several applications have been treated by Caputo 
himself up to nowadays, see e.g. [6l [7| and references therein. 



^According to Samko, Kilbas & Marichev^39 and Butzer & WestphalfS] the "regular- 
ized" fractional derivative was considered by Liouvillc himself (but then disregarded). 
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We point out that the Caputo fractional derivative satisfies the relevant 
property of being zero when applied to a constant, and, in general, to any 
power function of non-negative integer degree less than m , if its order jS is 
such that m — 1 < (3 < m . Furthermore, since 



we note that 



fit) = git) ^ fit) = git) + E c, t^-^ , (2.18) 

i=i 

whereas, using also (2.15) or (2.16), 

m 

,DUit) = tD^.git) ^ fit)=git) + Y.c,t"^-^. (2.19) 

In these formulae the coefficients cj are arbitrary constants. We also note 
the different behaviour of tD^ at the end points of the interval (m — 1, m) , 

hm ,Df/(t) = (t)-/(™-i)(0+), lim ,Df/(t) = /(-)(t). 

(2.20) 

The last limit can be formally obtained by recalling the formal representation 
of the m-th derivative of the Dirac function, 5'^^\t) = t~'"~^/r(— m) , t > , 
see [IE]. As a consequence of (2.20), with respect to the order, the Caputo 
derivative is an operator left-continuous at any positive integer. 

The Mellin transform 

Let 

M {fir); s} = f*is) = fir) r'-' dr, 71 < 5J (s) < 72 (2.21a) 

^0 

be the Mellin transform of a sufficiently well-behaved function /(r) , and let 
{f*is); r} = fir) = — / /* (s) ds (2.216) 
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be the inverse Melhn transforn|§, where r > , 7 = 3? (s) , 71 < 7 < 72 • 

We refer to speciahzed treatises and/or handbooks, see e.g. [12], [52] , 
|37j . for more details and tables on the Mellin transform. Here, for our 
convenience we recall the main rules that are useful to adapt the formulae 
from the handbooks and, meantime, are relevant in the following. 

Denoting by ^ the juxtaposition of a function /(r) with its Mellin trans- 
form f*{s) , the main rules are: 

fiar) ^ a-^ris), a>0, (2.22) 

r-f(r) ^ r(s + a), (2.23) 

00 

Hr) = [- fip) g{r/p) dp n h*is) = ris) g*is) . (2.25) 

The Mellin convolution formula (2.25) is useful in treating integrals of Fourier 
type for x = \x\> Q : 

1 

Ic{x) = — I f (k) cos (k x) dn , (2.26) 

TT JO 

1 r°° 

J,(x) = - / f{K) sin {Kx)dK, (2.27) 

71 Jo 

when the Mellin transform f*{s) of /(/t) is known. In fact we recognize 
that the integrals Ic{x) and Is{x) can be interpreted as Mellin convolutions 

^For the existence of the Mehin transform and the vahdity of the inversion formula we 
need to recall the following theorems, and adapted from Marichev's treatise [32], 
see THEOREMS 11, 12, at page 39. 

Let /(r) G L'^{e, -E),0<e<i?<oo,be continuous in the intervals (0, e] , [E, 00) , 
and let |/(r)| < M r'"'^ for < r < e , |/(r)| < M r-f^ for r > E , where M is a constant. 
Then for the existence of a strip in the s-plane in which f{r)r^^^ belongs to L^(0,oo) it 
is sufficient that 71 < 72 • When this condition holds, the Mellin transform /*(s) exists 
and is analytic in the vertical strip 71 < 7 = 5R(s) < 72 . 

(6b) jf jg piecewise differentiable, and /(r) r'''"^ g L'^(0,cx)) , then the formula (2.21b) 
holds true at all points where /(r) is continuous and the (complex) integral in it must be 
understood in the sense of the Cauchy principal value. 
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(2.25) between /(k) and the functions 5'c(k) , gs{n) , respectively, with r 
l/\x\ , p = K , where 



9c{k) := cos (- ) n sin ( ^ ) := g*{s), < 3?(s) < 1 



vr X K VK 



TT \X\ K, \K 



M 


r 


(1 






■ sin ( 




>• 




TT 












(1 




s) 


■ cos 1 


/vrs 




TT 


X 




It 



(2.28) 



g^{K) := sin (- ) n cos ( — ) := g*{s), < 3?(s) < 2 , 



(2.29) 

The Melhn transform pairs (2.28)-(2.29) have been adapted from the tables 
in [32] by using (2.22)-(2.24) and the duplication and reflection formulae for 
the Gamma function. Finally, the inverse Mellin transform representation 
(2.21b) provides the required integrals as 

1 1 ny+ioo / TV s\ 

Ic{x) = — — / r(s)r(l-s) sin — x^rfs, x>0, 0<7< 1, 

TTX iTXl Jf-ioo \ 2 / 

(2.30) 

1 1 /-l+ioo /vr s \ 

Isix) = — — r{s)T{l-s) COS — x^ds, x> 0, 0<7<2. 

TTX 2m J-y-ioo \ 2 / 

(2.31) 

3. Scaling and similarity properties of the 
Green function 

For the equation (1.1) we consider the Cauchy problem 

m(x,0) = (^(x) , xeiR, M(±oo,t) = 0, t>0, (3.1) 

where ip{x) G L'^{IR) is a sufficiently well-behaved function. If 1 < /3 < 2 we 
add to (3.1) the condition Mi(x, 0) = 0, where Ut{x,t) = ■^u{x,t). 

By solution of the Cauchy problem for the equation (1.1) we mean a 
function ^(x, t) which satisfies the conditions (3.1). By the Green function 
(or fundamental solution) of the Cauchy problem we mean the (generalized) 
function G^^(x, t) which, being the formal solution of (1.1) corresponding to 
ip{x) = 6{x) (the Dirac delta function), allows us to represent the solution 
of the Cauchy problem by the integral formula 

/ + 00 
G%{^,t)y,{x-Od^- (3.2) 
-co 
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It is straightforward to derive from (1.1) and (3.1) the Fourier-Laplace trans- 
form of the Green function by taking into account the Fourier transform for 
the i?ies2;-Fe//er space-fractional derivative, see (2.2), and the Laplace trans- 
form for the Caputo time-fractional derivative, see (2.12). We have (in an 
obvious notation): 

-rM ^A^^ = ^d^^ - ' (3-3) 

where 

:= |«re^(sig"'^)^^/2 =^;f(3;^ = ^-^(-«;). ^34) 
We therefore obtain 

By using the known scaling rules for the Fourier and Laplace transforms, 

f{ax) ^ a-^f{K/a), a>0, f {ht) ^ h'^ f{s/h) , b>0, (3.6) 

we infer directly from (3.5) (thus without inverting the two transforms) the 
following scaling property of the Green function, 

G%{ax , bt) = h-^G%{ax/h^ , t) , 7 = . (3.7) 

Gonsequently, introducing the similarity variable x/P , we can write 

GiA^, t) = t-^ K%{x/t^) , 7 = , (3.8) 

where the one- variable function K^^^ is to be determined as indicated below. 

Let us first invert the Laplace transform in (3.5). To this purpose we 
recall the Laplace transform pair, 

Ep{ct^) ^ 5ft(s)>|c|V^ (3.9) 

with cG(r, 0</9<2, where E/3 denotes the entire transcendental function, 
known as the Mittag-Leffler function of order (3 , defined in the complex plane 
by the power series 

00 ~n 
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Then, by comparing (3.5) with (3.9), we obtain the Fourier transform of the 
Green function as 



(3.11) 



For detailed information on the Mittag-Leffler-type functions and their 
Laplace transforms the reader may consult e.g. [11], [10], [39], [26], [20] , 
[21], [36]. We note for later use that the Mittag-Leffler function (3.10) ad- 
mits a Mellin transform type representation, originally due to Barnes [1] (see 
also [32] p. 118 (7.79)), as 



1 

27ri 



T{s)T{l-s) 



-z 



' ds . 



(3.12) 



where the integration is over a left-hand loop L_oo drawn round all the left- 
hand poles s = 0, —1, —2, ... of the integrand in a positive direction. Indeed, 
by the residue theorem it is not difficult to write the integral in (3.12) as the 
power series in (3.10). As a matter of fact, from (3.12) we can deduce the 
following Mellin transform pair (see also [32], p. 300) which will be used to 
get the Fourier antitransform of (3.11), 



M 



r(^)r(i-5) 



r>0, 0</5<2, 0<3fJ(s)<l. (3.13) 



In view of the symmetry properties of V'q(^) ? see (3.4), and of the Mittag- 



Leffier function i.e. Ep{z) = Ep{z) , z G(P , we have 



G%{K,t)=G%{-^,t)=G~'f,{-^,t). 
Furthermore, we easily recognize from (3.10)-(3.11), 

Gf^(0,t)=E^(0) = l, t>0. 



(3.14) 



(3.15) 



Provided that (^^^(x, t) does exist as inverse Fourier transform of (3.11), 
equations (3.14)-(3.15) ensure that ^^^(x,^) is real and normalized, i.e. 

The inversion of the Fourier transform (at most as a improper integral) 
requires (in view of the Riemann-Lebesgue lemma) that for t > , 



E, 



as \k\ 



oo . 



(3.16) 
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Taking into account the growth properties of the Mittag-Leffler function, see 
e.g. Gorenflo, Luchko and Rogosin [2D], this means that the LHS of (3.16) is 
required to be bounded for k E M and t > . It turns out that for this we 
must require 

\e\<2-f3. (3.17) 

This means that for 1 < /? < 2 the allowed {a , 9} region could be no longer 
the Feller- Takayasu diamond \6\ < min{a, 2 — a} ; in fact, in the cases 
0<a<l</3<2 with a + /3 > 2 and l<a</3<2we must cut-off the 
upper and lower corners of the diamond with the lines 6 = ±(2 — j3) . 

As far as the determination of the one- variable function K^p{x) is con- 
cerned, see (3.8), we note the symmetry relation 



As a consequence, we can restrict our attention to x > , and obtain 
1 



211 

rK. 



-IKX 



1 



E, 



73 



a,l3 



X = - 



,KA^) 



TT Jo 
1 



TT JO 



COS (kx) 5R 



sin (kx) $5 



KlAx) +sK^[x] 



i9n/2 



i97r/2^j 



a,l3\ 

dn , 
dn . 



(3.18) 

(3.19) 

(3.20) 
(3.21) 



From (3.19) we can easily obtain the value attained by K^^ p{x) in x = . To 
this purpose we extend the argument in |T7] writing 



KM 



TT 

COS 



|^"i?^(^-K"e^^^/2) dK 



2a 



1 

Tia Jo 



dr . 



(3.22) 



Thus, the last integral above is the Mellin transform of the Mittag-Leffler 
function at the point s = 1/a ; it turns out to be convergent under the 
conditions (3 = 1 with a > , and 0<l3<2,i3y^l with a > 1 . For the 



finite value of K^j^{x) at x 



we thus use (3.13) obtaining 



KM 



f 1 r(i/a)r(i - l/g) 

na r(l — P/a) 

— r 1/a cos — , 
ira \ 2a 



cos 



'On' 
2a 



1 < a < 2, /? ^ 1, 
< a <2, (3 = 1. 



(3.23) 



We note that JO) is non negative except for 1 < a < /3 < 2 . 
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4. Particular cases for the Green function 



In this Section we are going to consider the important particular cases of our 
space-time fractional diffusion equation, i.e. 

{a — 2 , (3 = 1} {standard diffusion), 

{0 < q; < 2 , ;9 = 1} {space fractional diffusion), 

{q; = 2,0</?<2,/?7^1} {time fractional diffusion), 

{0 < a = (3 < 2} {neutral fractional diffusion), 

It is well known that for the standard diffusion {a = 2 , (3 = 1} the Green 
function is the Gaussian pdf 

=r^/^--^exp[-xV(4i)] , -oo < x < +oo , t>0, (4.1) 
2^y^^ 

with similarity variable x/t^^"^ , that evolves in time with moments (of even 
order) 

li2n{t)-^ x^''GlAx,t)dx^^-f^f^ , n = 0,l,2,..., i>0. 

(4.2) 

The variance o"^ := ^2{t) = 2t is thus proportional to the first power of time, 
according to the Einstein diffusion law. 

In Fig. 2 we report the plots of the Gaussian. K2i{x) — 
l/(2y^) exp(— x^/4) in the interval — 5 < x < 5 by adopting for the or- 
dinate a hnear scale (left) and a logarithmic scale (right). We recognize that 
the logarithmic scale is to be preferred to point out the tails. At the end 
of this paper we shall exhibit a number of lin-log plots of the fundamental 
solution for different values of the parameters a, 6, {3, that the reader can 
compare with the corresponding plot of the Gaussian, and emphasize the 
different behaviour of the tails. 
The space-fractional diffusion 

Let us now consider {0 < a < 2 , /3 = 1} {space fractional diffusion including 
standard diffusion for a — 2). In this case, reducing the Mittag-Lcfilcr 
function in (3.11) to the exponential function, we recover the characteristic 
function of the class of Levy strictly stable densities according to the Feller 
parameterization, see (2.2), (2.3) and (3.4). In fact, denoting this class by 
{L^(x)} , we have 

l|(/^) =e-^'('^), and gJi(«:, i) = e -^^"('*) , (4.3) 



15 




-5 -4 -3 -2 -1 1 2 3 4 5 -5 ^-3 -2 -1 1 2 3 4 5 



Fig. 2 The Gaussian pdf. 



with ^^(k) given by (3.4) and a and 9 restricted in the Feller-Takayasu 
diamond, see (2.3). Then the Green function of the space-fractional diffusion 
equation can be interpreted as a Levy strictly stable pdf, evolving in time, 
according to 



Gii{x,t)^t-^^"Li{x/t^/''), -oo<x<+oo, t>0. 



(4.4) 



A stable pdf with extremal value of the skewness parameter is called 
extremal. One can prove that all the extremal stable pdfs^ with < a < 1 
are one-sided, the support being Mq if ^ = —a , and Mq if = +a . The 
one-sided stable pdf^s with support in Mq can be better characterized by 
their (spatial) Laplace transforms, which turn out to be 

roo a. 

L-"is):^ e~^^ L-"(x)dx^e~^ , 3? (s) > , < a < 1 . (4.5) 

^0 

The stable densities admit a representation in terms of elementary func- 
tions only in the following particular cases 
a — 2 , 9 — , Gauss : 



1 



20F 



-x74 



— OO < X < +00 



a— 1/2, 9 — —1/2, Levy-Smirnov : 

3.-3/2 



a — 1 



e * ^ ^ Ly^\x) 
9 — , Cauchy : 



2^ 



-\/{^x) 



a; > 0: 



L\{x) = - 



1 1 



TT -|- 1 ' 



-00 < X < +00 . 



(4.6) 



(4.7) 



(4.8) 
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We note that the case a = 1 can be easily treated also for 6^0 taking into 
account elementary properties of the Fourier transform; we have 
a = 1, < 1^1 < 1 : 



1 cos(^7r/2) 
n [x + sin(07r/2)]2 + [cos(^7r/2)]^ 



^i(^) = r , ^.^/zi,/oM2 , lL.rn^,n^^2 ' -oo < x < +00 ; (4.9) 



a = l, = ±1 : 

Lf\x) = S{x±l) , -cx)<x<+oo. (4.10) 

We note that in the above singular cases a = 1 , 6* = ±1 , the corresponding 
Riesz- Feller derivatives xDj.i reduce to ±xD, see (2.10), so our fractional 
diffusion equation (1.1) degenerates into kinematic {i.e. first-order) wave 
equations and the corresponding Green functions are Gf l{x,t) = 6{x ± t) , 
meaning pure drift. 

For < a < 2 the stable prf/'s exhibit fat tails in such a way that their 
absolute moment of order u is finite only if — 1 < i/ < a . In fact one can 
show that for non-Gaussian, not extremal, stable densities the asymptotic 
decay of the tails is 

0(|x|-("+^)) , x^±oo. (4.11) 



For the extremal densities with a ^ \ this is valid only for one tail, the 
other being of exponential order. For < a < 1 we have one-sided pdf^s: 
for 6 = —a the support is 1R~^ and the pdf tends exponentially to zero as 
X O''" ; for = +a the support is 1R~ and the pdf tends exponentially to 
zero as X ^ 0^ . For 1 < a < 2 the extremal pdf^s are two-sided and exhibit 
an exponential left tail (as x —>■ — oo) if ^ = +(2 — a) , or an exponential right 
tail (as X —>■ +oo) if 6 = — (2 — a) . Consequently, the Gaussian distribution 
is the unique stable distribution with finite variance. Furthermore, when 
< a < 1 , the first absolute moment is infinite so we should use the median 
instead of the non-existent expected value. However, there is a fundamental 
property shared by all the stable prf/'s that we like to point out: for any a 
the corresponding stable pdf is unimodal and indeed hell-shaped, i.e. its n-th 
derivative has exactly n zeros, see Gawronski |15j. 

A general representation of all stable pdf^s in terms of higher transcen- 
dental functions has been achieved only in 1986; it was Schneider [13], who 
first has proved that the stable pci/'s can be expressed by means of Fox H- 
functions, see also Uchaikin & Zolotarev 
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We note that already in 1952 Feller [13] had obtained, by inverting the 
Fourier transform of the characteristic function, representations of the stable 
pci/'s in terms of convergent and asymptotic power series. Feller's results 
have been revisited (and corrected) by Schneider and can be summarized as 
follows. We restrict our attention to x > , since the evaluations for a; < 
can be obtained using the symmetry relation L^{—x) = L~^{x) . 

The convergent expansions are 
< a < 1 , |6'| < a : 



1 < a < 2, 



1 -a^n^(l + n«) . 

> {-X j — ; sm 

TT X ~, n'. 



TTX 



n=l 

< 2 - a 

oo 
n=l 



{9 -a) 



a;>0; (4.12) 



T{l+n/a) 



sm 



nn 
2^' 



— a] 



X > 0. 



(4.13) 



From the series in (4.12) and the symmetry relation we note that the extremal 
stable distributions for < a < 1 are unilateral, that is vanishing for x > 
if 6' = a , vanishing for x < if ^ = —a . 

The asymptotic representations are given by 
0<q;<1, -a < 9 < a : 



Li{x) sin 



n=l 

0<q;<1, 9 = -a 



nir 
2^' 



— a 



x^O^ 



(4.14) 



X 



0+, = {[27r(l-«)]"V^/(^-")} 



1/2 



2- a 



(4.15) 



ai 



2(1 -a 

l<a<2, a-2<9 <2-a : 



6i = (1 - a) a"/(i-°) , ci 



a 



1 — a ' 



1 



r(l + na) 



-X 



sm 



TT X 



n=l 



nl 



nn 

T' 



a 



, X oo , (4.16) 



l<a<2, 9 = a-2 



LT^(x) ~ Aox^' e"^^"'' 



X 



oo , Aq 



2TT(a — 1) a 



"1/2 



(4.17) 
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2{a-l) ' ' ' a-1 

From (4.12), (4.13) and (4.14), (4.16), we thus note that for non extremal 
densities the roles of convergent and asymptotic power series are interchang- 
ing with respect to the cases < a < 1 and 1 < a < 2 . 

As a consequence of the convergence of the series in (4.12)-(4.13) and 
of the identity L^(— x) = L~^{x) . we recognize that the stable prf/'s with 
1 < a < 2 are entire functions whereas the stable prf/'s with < a < 1 have 
the form 

forx>0, , . 

\ (l/|x|)$2(|x|-") forx<0, ^ ' 

where $i(-2) and $2(-2) are entire functions. The case a = \ {\Q\ < 1) can be 
treated in the limit for a 1 of (4.12) and (4.13), where the corresponding 
series reduce to geometric series in 1/x and x , respectively, with a finite ra- 
dius of convergence. The corresponding stable pdf^s are no longer represented 
by entire functions, as can be noted directly from (4. 8)- (4. 9). 

From (4.12) -(4.13) a sort of reciprocity relationship between stable pdf^s 
with index a and 1/a can be derived as noted by Feller [H]. Assuming 
1/2 < a < 1 and x > , we obtain 

^L?/Jx-")=Lr(x), r = a(^ + l)-l. (4.19) 

A quick check shows that 6* falls within the prescribed range, |^*| < a, 
provided that \6\ < 2 — 1/a . 

At the end of the paper we shall exhibit some plots the fundamental 
solution of the space-fractional diffusion equation, 

<,(x,l)=<,(x)=L^(x) 0<a<2, (4.20) 

in the range |x| < 5 . 

The time- fractional diffusion 

Let us now consider the case {a = 2, < /3 < 2} {time-fractional diffusion 
including standard diffusion for (3 = 1) for which (3.11) reduces to 

Gf^(K,t) = (-K^t^) , keR, t>0. (4.21) 
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Following Mainardi, see e.g. [271 EH], the problem can be treated with the 
equivalent Laplace transform 

= ^//'-^e~l^l^''^', -oo<x<+oo, 3?(s)>0, (4.22) 

with solution 

G\p{x, t) = {\x\/tf^^^) , -oo < X < +00 , t > , (4.23) 

where Ma/2 denotes the so-called M function (of the Wright type) of order 
P/2 , whose general properties are briefly discussed below. 

The function M^{z) is defined for any order u G (0, 1) and V2; by 

MJz) = y , , ^ \ = -y ) '—-r(un) sm(-Kun) , (4.24) 

It is a special case of the Wright function on which the interested reader can 
inform himself in several books and articles, e.g. [TT], [2B],[TH1 HH]. It turns 
out that Miy{z) is an entire function of order p = 1/(1 — z/) , which provides 
a generalization of the Gaussian and of the Airy function. In fact we obtain 

Mi/2(^) = ^ exp (- z'/a) , My,{z) = 3'/' Ai {z/S'/') . (4.25) 

For our purposes (time-fractional diffusion equation) it is relevant to con- 
sider the function M^, (0 < z/ < 1) for positive argument. In the following 
we shall denote generically by r a positive variable (it can be t or the 
similarity variable |x|/t^/^) and by c a positive constant which plays the role 
of a scaling parameter. We then point out the related Laplace transform 
pairs, see e.g. Mainardi and Gorenfio, Luchko & Mainardi |18j . 

' M^{r/c) <^ cE^{-cs) , 3fJ(s)>0, (a) 

^ {cr-") ^ s''-^ expi-cs") , 3fJ (s) > , (b) (4.26) 
^ (cr-'^) h exp(-cs^) , 3? (s) > . (c) 

We note that, thanks to Bernstein theorem, Eq. (4.26-a) proves the pos- 
itivity of M^{r) for r > since for < u < 1 the function E^, is known 
to be completely monotonic on the negative real axis. Furthermore, in the 
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(singular) limit for v — > 1~ , Eq. (4.26-a) shows that My{r) tends to the 
Dirac function b{r — 1) . Mainardi, see e.g. EB], has provided plots for 
r > 0, which show that M^{r) is indeed positive, and is monotonically de- 
creasing for < u < 1/2, while for l/2<z/<lit exhibits a maximum 
whose position tends to r = l as consistently with the limiting 

case Mi(r) = 6{r — 1) . We recognize that the result in (4.23) obtained as 
Laplace inversion of (4.21) is included in the Laplace transform pair (4.26-b) 
by taking u = (3/2, c = \x\ and r = t . 

For the function My{r) it is worth to recall the asymptotic representation 



M^(r) ~ Ao^^-^/^exp (-F) , r ^ oo , 
= ^- , Y =il-u) (u" r)i/(i-'^) 



(4.27) 



and the integral 



M,ir)dr= S > -I . (4.28) 



As a consequence of Eqs. (4.23) and (4.27) G2^ij{x,t) with < /? < 2 can 
be interpreted as a symmetric spatial pdf evolving in time, with a stretched 
exponential decay. More precisely, we have 

G%{x, 1) = ^ M^/2(|x|) ~ Ax" e-'^' , x^+oo, (4.29) 



A = {2n{2 - (3) 2/3/(2-/5)^(2-2/3)/{2-/3) j 



2/5 2 , ^^2-2/(2-/3)^/3/(2-/3) 



-1/2 



(4.30) 



2(2-/5)' ' ' 2-/3 



Furthermore, using (4.28), the moments (of even order) of G2i3{x,t) are 



/i2n(t) := / G%{x, t) dx = ) ^ t^" , n = 0, 1, 2, . . . , t > . 

(4.31) 

It is interesting to compare this expression with the analogue (4.2) for the 
Gaussian. In particular, the variance := fi2 = 2t'^/r(/5 + l) is now propor- 
tional to the /3-th power of time, consistent with anomalous slow diffusion 
for < (3 < 1 and anomalous fast diffusion for 1 < (3 < 2 . 
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In the slow diffusion case {0 < (3 < 1) the pdf attains its maximum value 
at X = (where the first derivative is discontinuous) and exhibits exponential 
tails fatter than the Gaussian; in the fast diffusion case {1 < (3 < 2) the pdf 
attains two symmetric maxima that move apart from the origin with time 
and exhibits exponential tails thinner than the Gaussian. 

We also outline some other interesting properties of the M function, which 
put it in relationship with the class of extremal stable distributions. In fact, 
we note that the Laplace transform pair (4.26-c), compared with (4.5), allows 
us to relate the one-sided, extremal. Levy stable pdf of index z/ (0 < z/ < 1) 
with the M function of order u ; we have 

jj^L-^{^)=^^M'^{^) ^ 0<.<1, OO, r>0. (4.32) 

Incidentally, the above relation turns out to provide an alternative proof of 
the non-negativity of the M function for positive argument. Furthermore, 
putting in Eqs (4.28) and (4.32) = a , c = t > and r = x > , we get 
the Laplace transform pair related to the (spatial one-sided) extremal stable 
pdf (evolving in time) of index a (0 < a < 1), which represents the Green 
function of the space-fractional diffusion equation with 6 = —a , see Eqs. 
(4.4)-(4.5), and henceforth, 

e-*^" A t-'/"L-"{x/t'/") = G^;i{x, t) = c^^ M^{t/xn , (4.33) 

with 0<a<l, x>0, t>0. 

Then, in virtue of the Feller reciprocity relation (4.19) and using Eqs. 
(4.4) and (4.23), after some manipulations we can recover the noteworthy 
result 

Glp{x,t) = ^Glf^~^{x,t) , 1<(3<2, x>0,t>0. (4.34) 

This result (already noted by Mainardi & Tomirotti [3Tj by comparing the 
corresponding series expansions, see (4.13) and (4.24)) states that the two 
symmetric branches of the Green function of the time-fractional diffusion 
equation of order 1 < j3 < 2 are proportional to the corresponding (not 
symmetric) branches of the Green function of the space- fractional diffusion 
equation of order l<a = 2/P<2 with skewness 9 = ±(2 — /?/2), namely 
the exponential queues of the two extremal stable distributions of index 2/(3 . 
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In the limit /3 = 2 we recover the fundamental solution of the D'Alembert 
wave equation, i.e. 



Gl,(,, t) = + " = + , (4,35) 

where — oo < x < +00 , t > . 

At the end of the paper we shall exhibit some plots the fundamental 
solution of the time-fractional diffusion equation, 

G%{x,l) = Klp{x) = ]^Mp,^{x), 0</3<2, (4.36) 

in the range \x\ < 5 . 

The neutral-fractional diffusion 

Let us now consider {0 < a = /5 < 2} {neutral fractional diffusion)., which 
includes the Cauchy diffusion for a = (3 = 1 (0 = 0) and the limiting case of 
wave propagation for a = /5 = 2. 

For 1 < a < 2 and = the Green function has been derived in explicit 
form by Gorenfio, Iskenderov and Luchko [T7] by using the Mellin-Barnes 
integral representation. We shall adopt this representation later. Now, we 
consider it conceptually more economical to remain (as long as possible) in 
the kingdom of Fourier-Laplace transforms and we point out the following 
Fourier transform pair related to the Mittag-LefHer function of our interest: 

/ I i«x .7^ 1 1x1"""^ sin (Q;7r/2) ^ „ 



(4.37) 

This pair can be verified as an exercise in complex analysis following the 
method illustrated by Gorenfio and Mainardi, see e.g. \2lf[ As far as we 
know, this case of fractional diffusion seems not be treated in its greatest 
generality in the literature. Now, taking into account elementary properties 



In [21] the reader can find the proof of tlic Laplace transform pair 

^ ' TT 1 + 2i'^ cos (i/tt) + i2'> ' ' ' V / 



23 



of the Fourier transform and (4.37), the reduced Green function in (3.19) 
reads for a = /3 and x > as 



(x) = N^{x) = ^ — — , 0<a<2. (4.38) 

"'"^ ^ ^ TT l + 2x"cos[f(a-^)] +x2" ' ^ ' 

This solution, that can be extended to negative x by setting A^^(— x) = 
N~^{x) , is evidently not negative in all of M, so it can be interpreted as a 
probability density. In other words, N^{x) may be considered the fractional 
generalization with skewness of the well-known Cauchy density (4.8). In the 
limiting case a ^ 2^ (with 6 = 0) the density tends to the combination 
[6{x — 1) + 6{x + l)]/2 , so we recover the Green function of the D'Alembert 
wave equation quoted in (4.35). 



5. Composition rule for the Green function 
with < (3 <1 

We now present a composition rule which allows us to express the general 
Green function of the space-time fractional diffusion equation (with the re- 
striction < /5 < 1) as an integral involving the two Green functions corre- 
sponding to space-fractional and time-fractional diffusion equations. To this 
purpose we note that the Fourier Laplace transform of the Green function 
(3.5) can be re-written in integral form as in 



(5.1) 



(k si - ^'^ ^ - s^^-l n p-^i-^^ + ^al'^)] rill 



oo 

e' 



(^s^-ig-^^^) du. 
In view of Eqs. (4.3) and (4.26) we can interpret the above formula as 

roc - — . 

G%{k, s) = 2j^ Gl,{^, u) GO 2^(«, s) du . (5.2) 
Then, by inversion, we obtain the required composition rule 

roc 

Glpix, t) = 2 Gl,{x, u) Gl,p{u, t) du . (5.3) 
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Note the presence of (^2,2/3 instead of G^,/? • Hence Eq. (5.3) is a formula 
separating variables. It states that the Green function for the space-time- 
fractional diffusion equation of order {a, (3} , with < a < 2 and < /3 < 1 , 
can be expressed in terms of the Green function for the space-fractional 
diffusion equation of order a and the Green function for the time-fractional 
diffusion equation of order 2(3 . 

We now present alternative, equivalent forms of the composition rule, 
that directly involve the functions L and M , limiting ourselves to x > . 
Because of Eqs. (4.4) and (4.23), we can write Eq. (5.3) as 

GU^, t) = t-^ ru-'/- Ll (^) M, (^) du . (5.4) 



Remark 5.1. The formulae (4.4) and (4.23) corresponding to the particular 
cases {a, 1} {space-fractional diffusion) and {2,/5} {time-fractional diffusion) 
are recovered from (5.4) as follows: 



00 



(5.6) 

Eq. (5.6) is of high interest for the theory of the M functions since it is a 
sort of (integral) duplication formula with respect to the order; it is worth 
to note the presence of the Gaussian with a spreading variance in the kernel 
of the integral. 

Taking into account the relation (4.32), the function M in (5.4) can be 
expressed in terms of an L function. We have 



Gl,{x.t) 



u 



X 



u 



1/a 



du. (5.7) 



Putting y = u we derive from (5.7) the composition rule in the form 
recently obtained by Uchaikin, see |l6l HT], (for t = 1) : 

G%{xA) = KAx) = Hy^'-Li L/{y)dy. (5.8) 



The composition rule, in its equivalent forms (5.3), (5.4), (5.8), shows 
that the Green function of the space-time fractional diffusion equation is non 
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negative for any x e IR and t e IRq , when {0<q;<2, < (3 < 1} . We 
also know that this is the case when {a = 2 , < f3 < 2} {time-fractional 
diffusion) and when {0 < a = f3 < 2} {neutral-fractional diffusion). In 
the next Section, by using the Melhn convolution, we shall derive a new 
composition rule which ensures the non- negativity of K^i^{x) in all of M 
for {0 < q; < 2} n {0 < /9/q; < 1} . Then we can conclude that the Green 
function can be surely interpreted as a spatial pdf, evolving in time, in the 
ranges {0 < a < 2} n {0 < < 1} and {0 < /3 < a < 2} . 

6. Mellin-Barnes integral representations for 
the Green function 



Let us now consider the Fourier representation of the general Green function 
K^p{x) (restricting to x > 0) as stated at the end of the Sect 3, see Eqs 
(3.18)-(3.21). We intend to use the Mellin convolution to invert the relevant 
Fourier transforms according to the scheme (2.26)-(2-31). To this purpose 
we also need the Mellin transform pair deduced from (3.12)-(3.13) 

where k>0, |0|<2-/3, 0< U{s) < a . 

Using Eq. (6.1) in Eqs (2.28) and (2.29), Eq. (3.12) yields 

<^(x) = ,K%{x) + sK%(x) , x>0, (6.2) 

with 

, , 1 1 r+"~r(f)r(i-f)r(i-s) /irs\ _ /fas' 



' ' irax 2m Jj-ioo r(l - ^s) \2 J \2a J 

where < 7 < min{Q;, 1} , and |^| < 2-/3. We thus obtain from Eqs 
(6.2)-(6.3) 



. N _ 1 1 r+^°°r(^)r(i-^)r(i-s) . 
'^'^^'^^ " ^2^iJ,-ioo r(i - ^s) 



x^ ds. 
(6.4) 
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By setting 



P 



a-e 

2a 



(6.5) 



and using the reflection formula for the gamma function, we finally obtain 

1 1 n+i<^ r(^)r(i - ^)r(i - 



ax 2Tri 



■y—ioo 



r{l- ^s)T{ps)T{l-ps: 



ds . 



(6.6) 



By using a standard notation for the "gamma fraction", see e.g. [32], p. 65 
(4.13), we can re-write (6.6) as 



1 1 /•7+ioo 
J^—ioo 



''^ ax 2m 



- 1 - - 1 



1 



ps , 1 — ps 



x^ ds . 



(6.6') 



The formulae (6.4) and (6.6) are valid for < 7 < min{a, 1} with |^| < 2-/3 , 
and provide equivalent integral representations of the Green function for the 
general space-time fractional diffusion equation. In the following we shall use 
(6.4) or (6.6) [(6.6')] according to our convenience. 

The integral at the RHS of Eq. (6.6) is a particular Mellin-Barnes inte- 
gral according to a usual terminology, see e.g. [11], Vol. 1, Ch. 1, §1.19, 
pp. 49-50. The interested reader can find in pJ!] the discussion on the gen- 
eral conditions of convergence for the typical Mellin-Barnes integral, formerly 
given by Dixon & Ferrar [5], and based on the asymptotic representation of 
the gamma function. By using these results we can verify once again the 
convergence condition (3.17), namely |^| <2 — (3 . 

In the particular cases of standard diffusion, space-fractional diffusion, 
time-fractional diffusion and neutral-fractional diffusion, the Mellin-Barnes 
integrals in Eq. (6.6) simplify as follows. 

For the standard diffusion , namely for {a = 2 , 9 = , P = 1} , we have 
p = 1/2 and 



Kl{x) = Ll 



X 



1 1 
2x27ri 



7+ioo r(i-s) 



7—100 



r(i - s/2) 



X 



ds, < 7 < 1 . (6.7) 



^The names refer to the two authors, who in the first f 9f O's developed the theory of 
these integrals using them for a complete integration of the hypergeometric differential 
equation. However, as pointed out by Tricomi in [TT (Vol. 1, Ch. 1, §1.19, p. 49), these 
integrals were first used by S. Pincherle in 1888. For a revisited analysis of the pioneering 
work of Pincherle (1853-1936, Professor of Mathematics at the University of Bologna from 
1880 to 1928) we refer to the recent paper by Mainardi and Pagnini [3D]. 
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For the space-fractional diffusion ({0 < o; < 2, |^| < min{Q;, 2 — a}, 
/3 = 1}) we distinguish a number of cases. 
For {0 < q; < 1 , 1^1 < q; , /9 = 1} we have < p < 1 and 

Kl,{x) = Llix) = —— ) ) i-x'ds, 0<7<«. 

ax 2m J-y-ioo L [p s) L [1 — ps) 

(6.8) 

Let us consider the extremal cases. For = —a , we have p = 1 so r(l — 
s) cancels with r(l — ps) . Then the "gamma fraction" in (6.8) reduces to 
r(s/Q;)/r(s) . For ^ = a , we have p = so K^^^ix) = for x > . 
For {a = 1 , 1^1 < 1 , /3 = 1} we have < p < 1 and 

Let us consider the extremal cases. For ^ = — 1 , we have p = 1 so r(l — 
s) cancels with r(l — ps) and T{s) cancels with r(ps) . Then the "gamma 
fraction" in (6.9) reduces to 1 to yield the Mellin representation of the Dirac 
function S{x — 1) . For ^ = 1 , we have p = so Kl ^^{x) = for a; > . 
For {1 < a < 2 , \9\ < 2 - a , P ^ 1} we have 0<{a-l)/a<p<l/a<l 
and 

KM--Ll{x) = —— ) ' ' \ '-x'ds, 0<7<1. 

ax 2m J-t-ioo r(ps) r(l — ps) 

(6.10) 

Let us consider the extremal cases. For 9 — —{2— a) , we have p = {a—l)/a > 
so (6.10) is still valid. For 6 = 2 — a we have p = l/a so r{s/a) cancels with 
r(ps) , and the "gamma fraction" in (6.10) reduces to r(l — s)/r(l — s/a) . 

For the time-fractional diffusion ({a = 2, 9 = 0, Q<P<2, 1}) we 
put a = 2 and p = 1/2 in (6.6) and obtain 

(6.11) 

For the neutral-fractional diffusion ({0 < a = (3 < 2}) we put a = (3 in 
(6.6) and obtain 
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In the limiting case a = f3 = 2{6 = 0,p = 1/2) we recover from (6.11) 
or (6.12) the Melhn-Barnes representation of S{x — l)/2, consistently with 
the D'Alembert wave equation, see Eq. (4.35). We note that the Mellin 
transform pair (6.6) allows us to compute the value of any convergent moment 
of the Green function K^p{x) . In fact, recalling the basic formulae (2.21a)- 
(2.21b) for the Mellin transformation (2.21a)-(2.21b), from (6.6) we write the 
Mellin transform of xK^^{x) as 

KU-).U.-' r(-v«)r(i + .Va)r(i + .) 



"'^^ ' a r(l + /5s/a)r(-ps)r(l + ps) ' 

(6.13) 

-min{a, 1} < 3?(s) < . 

In order to include s = in the convergence strip (so, in particular, the 
integral of K^j^{x) in Mq can be evaluated) we properly use in (6.13) the 
known property r(l + z) = zr{z) to obtain 



r+oo 

Klp{x)x-'^dx = p 



r(l-s/a)r(l + s/a)r(l + s) 



T{l-ps)T{l + ps)T{l + l3s/a) ' 

(6.14) 

— minja, 1} < 9?(s) < a . 

In particular, setting s = we find K^^ p{x) dx = p (p = 1/2 if 6' = 0). 
Eq. (6.14) is consistent with a similar expression given by Uchaikin [17] . We 
note that it is strictly valid as soon as cancellations in the "gamma fraction" 
at the RHS cannot be done. Then Eq. (6.14) allows us to evaluate (in 
Mq) the (absolute) moments of the fundamental solution of order S such 
that — min{a, 1} < 6 < a . In other words, it states that, as x —>■ +oo, 
K^ij{x) = O (^x"*^""*"^-*^ . When cancellations occur in the "gamma fraction" 
the range of 6 may change. In this respect an interesting case is {a = 2 , 6 = 
, < (3 < 2} {time-fractional diffusion including standard diffusion), where 
Eq. (6.14) reduces to 

/■+°° 1 m -I- 

I KU.).U.= -^AJ^^, s,,,>_i. (6.15) 

This result is consistent with the existence of all moments of order 6 > —1 
for the corresponding Green function, see (4.28). 
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Furthermore, by using the Melhn convolution formula (2.25), and the 
Mellin-Barnes representations valid in the particular cases, namely (6.8), 
(6.10) for the space-fraction diffusion, (6.11) for the time-fraction diffusion 
and (6.12) for the neutral-fraction diffusion, Eq. (6.6) can help us to inves- 
tigate the non-negativity of the Green function ^{x) . In fact, it is not 
difficult to recognize that, assuming x > and < a < 2, 

a /"[r-^ (r )] Li (x/O f , < /3 < 1 , 

J^Mp/^iOKix/Oj, 0<(3/a<l. 

The limiting cases jS = 1, a = (3 and a = 2 can be included if one takes 
into account, in a proper way, the resulting expressions by the generalized 
function 6{x — 1) . We note that the first Mellin convolution in (6.16) is 
consistent with the composition rules obtained in Sect 5 by using the Fourier 
transforn|§. 

The formulae in (6.16), by involving two non negative functions (related 
to the probability densities L^x), M/3/2(x) or N^{x)), allow us to state the 
probability interpretation of the Green function in the ranges 

{0 < a < 2} n {0 < /? < 1} 

and 

{0 < a < 2} n {0 < < 1} . 

Therefore, the probability interpretation holds true for any a G (0, 2] if 
< /3 < 1 , whereas, if 1 < /3 < 2 , only for 1 < /3 < a < 2 . The cases 
excluded from the probability interpretation turn out to be 

{0 < a < /?} n {1 < /3 < 2}. 



^To this purpose, let us consider the composition rule in form (5.4), 

Putting ^ ~ u^/" and taking into account the scaling property G^^(x,i) = 
f-P/a ^ (x/i*^/") , we rewrite the composition rule (for t = 1) in the form of the first 
Eq. in (6!l6). 
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7. Computational representations for the 
Green function 

Those who are acquainted with Fox H functions can recognize in (6.6) the 
representation of a certain function of this clasi^ on which the interested 
reader can inform himself in several books, e.g. [33], [H], [22], [SZ], [SH], [2S]- 
Unfortunately, as far as we know, computing routines for this general class 
of special functions are not yet available. Here, following and completing the 
approach adopted by Gorenflo, Iskenderov & Luchko [17], we intend to com- 
pute the Green function in any space domain providing for it a computational 
representation (usually to be obtained by matching two distinct expressions). 

To this purpose we distinguish three cases depending on the order rela- 
tionship between a and (3. According to their increasing difficulty we consider 

(i) a = 13, (ii) a</3, (iii) a > {3 . (7.1) 

For the case (i) {a = (3} the integral representation (6.4) simplifies into 

x-'^ds. (7.2) 

In this case the contour of integration can be transformed to the loop L_oo 
starting and ending at infinity and encircling all the poles s„ = —an, n = 
0, 1, 2, . . . of the function T{s/a) for < a; < 1 and to the loop L^^o starting 
and ending at infinity and encircling all the poles Sn = {l+n)a, n = 0, 1, 2, . . . 
of the function r(l — s/a) for x > 1 . Applying the residue theorem we arrive 
at the series representations 

(-x°)", 0<a;<l; (7.3) 

(-a;-°)", l<a;<oo. (7.4) 

So we have two different representations by power series: the one in positive 
powers, convergent for < x < 1, the other in negative powers, convergent 

""^^ ADDED NOTE (2007) As a complement of the present analysis the reader can consult 
the paper by F. Mainardi, G. Pagnini and R.K. Saxena: Fox H functions in fractional 
diffusion, published in Journal of Computational and Applied Mathematics 178 (2005) 
321-331, where the Authors have provided the representation of the fundamental solutions 
of the fractional diffusion equations here treated in terms of Fox H functions 
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Ka,ai^) = — II sin 
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for X > 1 . Remarkably, in this special case, we can obtain the fundamental 
solution in all of iR in a closed form, namely expressed in term of elementary 
functions. In fact, following ^l7\, we use the formula 



oo 

r" sin(ria) 

n=l 
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re 



r sm a 



\n=l 



re' 



with a E M, |r| < 1 , and we get 



KM 



1 



x°' ^ sin[|(a — 



TT 1 + 2a;"cos[f (a 



+ X 



2a 



1 — 2r COS a + r^ 



< X < oo . 



(7.5) 



(7.6) 



This result is in agreement with (4.38) and, by taking into account the 
symmetry relation (3.18), can be properly extended to — oo < x < 0. We 
also note that from (7.6) we recover the expressions of the particular cases 
{a = P = 1 , \9\ < 1} given in Eq. (4.9), namely the stable density L^{x), 
and {l<a = l3<2,9 = 0} given by Gorenflo, Iskenderov and Luchko, see 
Eq. (20) in [r7|. In the limit as x — , we get 



lunKjx) 



lim 



X 



a-1 



sm 



TC 



2^" 




< a < 1 , 

a = l, (7.7) 

1 < a < 2 , 



which is consistent with (3.23). 

For the case (ii) {a < (3} the contour of integration in (6.6) can be trans- 
formed to the loop L_oo starting and ending at infinity and encircling all 
the poles Sn = —an, n G INq = {0, 1,2, . . .} of the function T{s/a) . The 
residue theorem for simple poles gives us the following representation by a 
convergent series in negative powers of x° , 



1 ^ r(l + an) . 
— > — sm 

7CX r(l + (3n) 



nil 

T' 



a] 



—X 



< X < oo. (7.8) 



We note that for f3 = 1 we recover the series at the RHS of (4.12) giving the 
representation of the stable density /^^(x) with < a < 1 . 

For the case (iii) {a > (3} the situation is more complicated because we have 
to transform the contour of integration in (6.6) to the loop Lj^^ encircling 
all the poles is„ = a(l + n), n E INq and = I + m, m E INq of the 
functions r(l — -) and r(l — s), respectively. 
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Wc note that for /9 = 1 or o; = 2 the poles i,s„ do not appear because 
r(l — s/a) cancels in (6.6) as seen in Eqs (6.7)-(6.11). 

Excluding the above cases, i.e. ii j3 ^ 1 and a 7^ 2 , we have to consider 
the possibility of double poles occurring when —a{k + l) = — (m + 1) , namely 
when a — {m + l)/{k + l) , m,k E Nq. 

When the poles are all simple the residue theorem gives us the following 
representation by two convergent series in positive powers of x°' and x, 



1 °° rfi 

k=0 



ak) 



nx " r(i 



Pk) 



sm 



a 



i-x-) 



(7.9) 
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TTX 



g r(i 



)r(i + 



sm 



kTT 

2^* 



— a) 



[-xf. 



,=0 ^!r(i - ^^k) 

We note that in the limiting cases /3 = 1 or a = 2 the first series has vanishing 
coefficients and the second series yields the expected results, namely the 
expansion in positive powers of x for the stable density L^^i.^) with a > 1, 
see (4.13), and for the density Mg/2(kl)/2, see (4.24) and (4.36), respectively. 

When the poles are not all simple, we need to separate the double poles 
from the simple ones. To this purpose we denote the szmj>/e poles of r(l— s/a) 
and r(l — s) by iSk = a{k + l),k E K , and = m + 1, m G M , and the 
double poles of r(l — -)r(l — s) by = (mo + 1)(1 + i),i G IVq where 



K = No\Ko; No = {0, 1, 2, ...}, Ko^{keN\k^ {ko+l)+{ko+l)i, i e TVq} , 

M = No\Mo; A/'o = {0, 1, 2, ...}, Mo = {m e iV I m = (mo+l) + (mo+l)i, i e No} 

Then, applying the residues theorem for simple and double poles we arrive 
at the following representation with three series (of which the first in positive 
powers of x", the second in positive powers of x, and the third of a peculiar 
character in that it contains positive powers of x that can be coupled with 
logx): 



1 ^ r{l -ak) . 

^ h r(i - m '''' 



^{9 
2 ^ 



a) 



-X 



a\k 



1 ^ r(i-^)r(i + ^) 



TTX 



m\T(l - ^m) 



sm 



ran 
2^* 



a 



(7.10) 

-x)^ + PU^). 
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Here P^jsi^) denotes the contribution from the double poles, which for a 
reads (see detailed derivation in Pagnini's thesis [55]): 



X sm 



E[(-i; 



X 



i=Q 



D 



(mp + 1)(1 + z)7r 
2a 



{9 -a) 



where 



Ri 

R2-- 



ip [(mo + 1)(1 + - logx, 
(mo + 1)(1 + i)7r 



a 

a) cot 



2a 



a] 



D = T[l- f3{h + 1)(1 + 0] r [(mo + 1)(1 + t)] 
and is the logarithmic derivative of the gamma function 

^(.)^-iogr(.)^^. 



(7.11) 



(7.12) 
(7.13) 

(7.14) 
(7.15) 

(7.16) 



Remark 6.1. For some rational values of the parameter (3, the relation 
1 — /3(/co + 1)(2^ + 1) = — /, Z = 0, 1, 2, . . . can take place for some values of 
the index i, z = 0, 1, 2 . . .. In this case, the term R2/D in (7.11) is an inde- 
terminate expression of the form — . Due to the formulae lims_»_i l/r(s) = 



0, / = 0, 1, 2 . . . , and lim. 



0, 1, 2, ... , we can 



rewrite Jx) in the form (see detailed derivation in Pagnini's thesis [35]): 



TT 



a 



^(_l)/3(feo+l)(l+i)+ 



i=0 



X sm 



(mo + l)(l + i)7r 
2a 



r[(mo + l)(l+2)] 

(7.17) 



— a 



which gives us a finite value. 

Remark 6.2. From the above analysis we recognize that all series, that are 
involved in the cases (i), (ii), (iii) . are vanishing if 6^ = a . Since this extremal 
value of 6 can be allowed only if < a < 1 , we can conclude that the Green 
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functions K"p{x) = for x > if a G (0, 1] and [3 G (0, 2) , subjected to the 
condition = a <2 — f3 . 

We must note that the representation of the Green function expressed 
for the cases (ii) and (iii) in terms of convergent expansions in negative and 
positive power series, respectively, are correspondingly not suitable to numer- 
ically evaluate the solution as soon as x is sufficiently small or sufficiently 
large. In other words, the convergent expansions cannot reproduce the ex- 
pected behaviour of their sum for x — > O"*" or a; — *^ oo when the corresponding 
series are in negative or positive powers of x , respectively. So, to complete 
our analysis we need to give the asymptotic representations of ii'^^(x) in the 
case (ii) {a < /?} as x — > 0+ and in the case (iii) {a > /5} as x — >• +oo . 
To this purpose we shall adapt to our analysis the results contained in the 
fundamental paper by Braaksma on the asymptotic expansions for certain 
Mellin-Barnes integrals. We shall limit ourselves to those cases (not already 
treated in Sect. 4) where the Green function is expected to be non-negative. 
As a consequence, for the case (ii) we consider < a < /9 < 1 , whereas for 
the case (iii) we consider < /5 < a < 2 with (3^1. 

For the case (ii) {{] < a < (3 < 1 , |^|<a} the asymptotic representation 
as X — > 0"^ turns out to be given by the convergent representation of the case 
(iii), namely by the two series in (7.9) or by three series in (7. 10)- (7. 17). We 
agree to write (for economy) 

In the extremal case 9 = —a only the first series in (7.18) [namely in (7.9) 
or (7.10)] survives, so we obtain the asymptotic representation 

<.W~-i: ^-0+. (7.19) 



X 



T{ak) r(l - (3k) 



We note that in the limiting case < a < (3 = 1 all the coefficients of the 
powers series in (7.19) are vanishing in view of the exponential decay just 
expected for the unilateral stable distribution L~"(x) as x — > 0+ . 

For the case (iii) {0 < (3 < a} we distinguish the sub-cases < a < 1 , 
where |^| < a , and 1 < a < 2 , where \6\ < 2 — a . In both sub-cases, when 9 
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does not assume extremal values, the asymptotic representation as x — ^ +00 
turns out to be given by the convergent representation of the case (ii), namely 



T{l + an) . 



sm 



—{e-a] 



-x"")", (7.20) 



As for the extremal cases, we have as follows: 

If < a < 1 , the representation (7.20) is still valid for ^ = — a , whereas 
is trivially zero for 9 — +a , as expected. If 1 < a < 2 , the representation 
(7.20) is still valid for ^ = 2 — a , whereas for ^ = a — 2 we have: 



Ax^e 



a „— bi'^ 



X 



+ 00 , 



(7.21) 



where 



A = {27r(a - I3)a^'^''-^^ 



/3/(«-/3)/?(«-2/3)/(a-/3)" 



-1/2 



(7.22) 



a = 



2(3 -a 



b={a- /?)a-"/("-^)/3^/("-^) , 



a 



2(q;-/3)' " '"^ "^'"^ ' ^ a-p 

Wc note that the above asymptotic representation (7.21)-(7.22) is still valid 
for /? = 1 , when it reduces to that for the density L"~^(a;) , see (4.17), and 
for q; = 2 (so 9 = 0), when it reduces to that for the density /^^^^(a;) = 
iM^/2(|a;|),see (4.29)-(4.30). 



8. Conclusive discussion and plots 

We can conclude with a discussion about some general features occurring 
in the Cauchy problem of our space-time fractional diffusion equation (1.1)- 
(1.2). A first general feature concerns the scaling property of the Green 
function which allows us to express it in terms of a function of a single vari- 
able, the reduced Green function K^ p{x), see (1.3) or (3.8). In this paper 
we have focused our attention to derive a computational form for ^{x) 
in all of M, taking into account the symmetry relation (3.18). The relevant 
particular cases of space-fractional ({0 < a < 2 , /? = 1}), time-fractional 
({a = 2 , < /3 < 2}) and neutral-fractional {{0 < a = (3 < 2}) diffusion 
have been summarized in Sect. 4 where the interpretation of the correspond- 
ing Green function as a probability density has been pointed out. 
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For computational aims the representation of ^{x) through the MeUin- 
Barnes integral, see (6.6), was found useful. Incidentally, this representation 
has enabled us to extend the probability interpretation of the reduced Green 
function to the ranges {0 < a < 2} fl {0 < < 1} and {1 < /5 < a < 2}. 

More precisely, to compute the function K^ p{x) we used the series expan- 
sions (7.8)-(7.17) and asymptotic expansions (7.18)-(7.22) which were derived 
from the representation (6.6). 

At the end of the section we shall exhibit some plots of the reduced 
Green function for some "characteristic" values of the parameters a, /3, and 
9. All the plots were drown by using the MATLAB system for the values 
of the independent variable x in the range |x| < 5. To give the reader a 
better impression about the behaviour of the tails the logarithmic scale was 
adopted. 

After long discussions we chose from the uncountable set of possible plots 
of the function K^ p{x) the plots having, in our opinion, the greatest interest, 
in the framework of the probability interpretation. 

I. For the Green function of the space-fractional diffusion equation {(3 — 1), 

we present plots for the cases: 
Fig. 3: a = 0.50 with a) ^ = and b) ^ = -0.50 ; 
Fig. 4 : a = 1 with a) ^ = and b) ^ = -0.99 ; 
Fig. 5 : a = 1.50 with a) ^ = and b) ^ = -0.50 . 

Note that for a = 1 we have chosen a nearly extremal case to show a density 
approaching to the delta function S{x — 1), see (4.10). 

II. The fundamental solution of the time-fractional diffusion equation {a — 2) 
was plotted in the following cases: 

Fig. 6: a) p = 0.25 and h) p = 0.50 ; 
Fig. 7 : a) p = 0.75 and b) /? = 1.25 ; 
Fig. 8: a) P= 1.50 and h) P = 1.75 . 

III. For the fundamental solution of the space-time-fractional diffusion equa- 
tion, we exhibit the following plots: 
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